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Abstract. We study a deterministic method for particle transport in tissue 
in selected medical applications. Generalized Fokker-Planck (GFP) theory |39| 
has been developed to improve the Fokker-Planck (FP) equation in cases where 
scattering is forward-peaked and where there is a sufficient amount of large- 
angle scattering. We compare grid-based numerical solutions to Fokker-Planck 
and Generalized Fokker-Planck (GFP) in realistic applications. Electron dose 
calculations in heterogeneous parts of the human body are performed. Ac- 
curate electron scattering cross sections are therefore included and their in- 
corporation in our model is extensively described. Moreover, we solve GFP 
approximations of the radiative transport equation to investigate reflectance 
and transmittance of light in tissue. All results are compared with either 
Monte Carlo or discrete-ordinates transport solutions. 



1. Introduction 

A difficult and important challenge in electron and photon transport is still the nu- 
merical solution of the Boltzmann transport equation (BTE) |5] . To contribute 
to this field of research, we study selected approximations of the transport equation 
for electrons and photons and present numerical results in different geometries. 

Nowadays cancer patients often undergo therapies with high energy ionizing 
radiation. In external radiotherapy photon beams dominate in clinical use whereas 
less patients receive electron therapy. Treatments are also performed by using 
heavy charged particles like ions or protons. These have higher costs for their 
particle accelerators but gain more and more importance due to first high intensity 
laser systems for protons [48j. 

To aid the recovery of patients it is important to deposit a sufficient amount 
of energy in the tumor. Simultaneously, the ambient healthy tissue should not be 
damaged. Therefore, the success of such radiation treatments strongly depends 
on the correct dose distribution. It is recommended that uncertainties in dose 
distributions should be less than 2% to get an overall desired accuracy of 3% in the 
delivered dose to a volume [H]. Additionally, thresholds have been developed to 
compare dose results computed by different algorithms. A suggested tolerance in 
homogeneous geometries is the 2% (relative pointwise difference) or 2 mm (absolute 
distance to agreement) criterion. However, in heterogeneities this limit increases to 
3% or 3 mm [54]. 

Up to now many clinical dose calculation algorithms rely on pencil beam mod- 
els. Originally developed for cosmic ray showers, Fermi ([T7], cited by |15]) and 
Eyges [TH] introduced a small-angle scattering theory which was afterwards applied 
to electrons. This theory was used by Hogstrom et al. to propose the pencil-beam 
model. Their algorithm includes experimental data, taken from dose measurements 
in a water phantom, to compute the central-axis dose [27]. Although it was a first 
clinically applicable model its accuracy deep in the irradiated material or in het- 
erogeneities is poor. This is basically due to the small-angle approximation in 
Fermi- Eyges theory, 9 « tan(0). It is true that single electron collisions show small 
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deviations. However, after multiple scattering events they accumulate to big angle 
changes in large penetration depths. Besides crude approximations like small devi- 
ations of particles throughout their whole path through the tissue [37], geometric 
structures transverse to the beam direction are assumed to be infinite. Although 
many improvements of Fermi-Eyges theory were performed, e.g., by including ad- 
ditional correction factors [33J EH HF" they still suffer from the small anlge and 
the homogeneity assumption. Comparisons with experimental data showed disad- 
vantages in inhomogeneous phantoms [41] . 

A statistical simulation method for radiation transport problems is the Monte 
Carlo method \4tj . It performs direct simulations of individual particle tracks which 
result from a random sequence of free flights and interaction events. In this way 
random histories are generated. If their number is large enough macroscopic quan- 
tities can be obtained by averging over the simulated histories. Monte Carlo tools 
model physical processes very precisely and can handle arbitrary geometries with- 
out losing accuracy. Although they rank among the most accurate methods for 
predicting absorbed dose distributions, their high computation times limit their 
use in clinics. Due to the increase in computing power and decrease in hardware 
costs Monte Carlo techniques have recently become a growing field in radiotherapy 
|51j . Not only general-purpose Monte Carlo codes are now publicly available but 
also commercial Monte Carlo treatment planning systems [501 HI]. However, they 
have not yet gained widespread clinical use. 

A different approach in the solution of radiation transport problems are de- 
terministic calculations solving the linear BTE. In principle, its solution will give 
very accurate dose distributions comparable to Monte Carlo simulations. The BTE 
can be analytically solved only in very simplified geometries, which is insufficient 
for clinical applications. Additionally, numerical solutions to the BTE require de- 
terministic methods coping with a six-dimensional phase space. Because of this 
and their simpler implementation, Monte Carlo techniques have so far prevailed 
in the medical physics community. Nevertheless, Borgers [8] argued that on cer- 
tain accuracy conditions deterministic methods could compete with Monte Carlo 
calculations. 

In this paper, we describe electron and photon transport in media by solving the 
Generalized Fokker-Planck (GFP) approximation of the linear BTE [SET. For elec- 
tron transport, this is an extension of the Fokker-Planck (FP) model, formulated 
in |25j . to higher-order approximations in angle scattering. It should be stressed 
that this approach brings along several advantages: The BTE does not require any 
assumption on the geometry so that arbitrary heterogeneities are possible. Further- 
more, we benefit from mathematical and physical approximation ansatze because 
they can be directly included in the differential equation. This avoids heuristic 
assumptions. In contrast to Monte Carlo simulation, deterministic solutions do not 
suffer from statistical noise and their resolution docs not depend on the number of 
particles traversing a certain region. Moreover, the treatment planning problem can 
be formulated as a PDE-constrained optimization problem. This structure can be 
used to obtain additional information to speed-up the optimization [T3 [20] • This 
is hard to achieve by Monte Carlo techniques. 

Previous studies in deterministic methods for radiotherapy primarily concen- 
trated on a combination of rigorous analytic solutions and laboratory measurements 
(see pencil beam models above). Without explicitly using experimental data in their 
model, Huizenga and Storchi [28j presented the phase space evolution (PSE) model 
for electrons and subsequently applied it to multilayered geometries [42] . Various 
improvements and extension to 3D beam dose calculations were performed [311132]. 
However, 3D dose calculations showed disadvantages in computation times [55] , 
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Moreover, PSE models use first-order discretizations in space which cannot com- 
pete with MC techniques [8]. By contrast, our access focuses on the continuous 
model of the Boltzmann equation discretized with high-order schemes. 

First studies for deterministic dose calculations in charged particle transport 
came from well-known procedures for numerical solutions to the transport equa- 
tion of neutrons or photons. Several approximative models to the BTE have been 
developed. Each of these methods has its advantages and drawbacks j9j. Multi- 
group methods are sometimes used to discretize the energy domain [13] . This leads 
to a number of monoenergetic equations to be separately discretized in space and 
angle. Whereas discretizations in space are usually done by finite difference and 
finite element (FE) methods, the remaining angle domain is discretized by discrete 
ordinates [6] [15] . Such an approach is implemented in the solver Attila [23] . Re- 
cently, 3D dose calculations for real clinical test cases were performed with Attila 
|53| . Results with similar accuracy as Monte Carlo calculations were achieved in 
promising computation times. 

Less frequent are angular FE approaches [TP seeking to reduce ray effects. Tervo 
et al. extended FE discretizations to all variables in spatial, angular and energy 
domains |52j so that no group cross sections were needed. Moreover, a detailed 
description of three coupled BTEs with FE discretizations of all variables was pro- 
posed in [7]. 

In this paper, the BTE is approximated by the continuous slowing down (CSD) 
method Scattering processes are modelled by Generalized Fokker-Planck ap- 
proximations [39] and compared with classical Fokker-Planck [l^ solutions. Leakeas 
and Larsen [391 showed that scattering kernels with a sufficient amount of large- 
angle scattering yield inaccurate FP results. As an extension of the work in S9J, we 
perform deterministic GFP simulations and investigate their behaviour in real ap- 
plications. It is well-known that electron scattering is dominantly forward-peaked. 
Hence, many electron transport simulations used the classical FP approximation 
|25| , I18| , [T4]. However, up to now no comparisons with GFP dose computations have 
been done including realistic physical scattering cross sections. In our case the lat- 
ter are extracted from ICRU libraries f25|. We describe in detail how transport 
coefficients in the BTE can be computed from these scattering cross sections and 
compute GFP electron dose profiles in inhomogeneous geometries. 

Even more challenging for GFP theory are scattering kernels including large 
angle-scattering. We therefore investigate transport of photons in tissues with 
forward-peaked and large-angle scattering. Using test cases from [24 the radiative 
transport equation for GFP approximations up to order five is solved to determine 
reflectance and transmittance of light in tissue. 

In the remaining paper you will find the following structure: A short discussion 
of basic electron interactions with matter is given at the beginning. We describe 
a model for electron transport and review crucial steps of the GFP theory. In 
addition, we derive transport coefficients for the GFP equations from databases 
for electron scattering cross sections. In Section 3, a deterministic model for light 
propagation together with different scattering kernels is introduced. Discretization 
methods used in our GFP equations are studied in Section 4. In Section 5, we com- 
pute FP, GFP and discrete-ordinates results for transmittance and reflectance of 
light by a slab. Moreover, we numerically compare FP, GFP and Monte Carlo solu- 
tions for 5 and 10 MeV electrons in homogeneous and heterogenous slab geometries. 
Section 6 gives the conclusions and outlooks. Appendix A contains explicit formu- 
lae for high-order polynomial operators. In Appendix B, we present the equations 
to be solved for the GFP coefficients. 
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2. Deterministic Model for Electron Transport 

2.1. Physical Interactions. Electron beams are nowadays a widely spread tool in 
cancer therapy. Typical electron beams, provided by high energy linear accelerators, 
range from 1 to 25 MeV. During irradiation of human tissue electrons interact with 
matter through several competing mechanisms: 

1. Elastic Scattering: This is usually a non-radiative interaction between electrons 
and the atomic shell. Projectiles experience small deflections and lose little 
energy. High energy electrons can also penetrate through atomic shells and are 
afterwards scattered at the bare nucleus without any energy loss. With kinetic 
energies above 1 keV elastic scattering in water dominantly occurs in the forward 
direction ^3S]. 

2. Soft Inelastic e^-e^ Scattering: Electrons interact with other electrons of the 
outer atomic shell which usually leads to excitation or ionisation of the target 
particle. Here binding energies are only a few eV so that projectile electrons 
transfer little energy and are hardly deflected. 

3. Hard Inelastic e~ -e~ Scattering: These collisions are determined by large trans- 
fer energies to the target electron. What 'large' exactly means is specifled in MC 
codes by cutoff energies. In PENELOPE [57|, for example, the default value of 
this simulation parameter is set to 1 % of the maximum energy of all particles. 
As a consequence, the target electrons are ejected with larger scattering angles 
and higher kinetic energies (delta rays) . They act as an additional source in the 
transport equation. 

4. Bremsstrahlung: Caused by the electrostatic fleld of atoms, electrons are accel- 
erated and hence emit bremsstrahlung photons. However, for energies below 
1 MeV this phenomenon can be neglected. Bremsstrahlung photons are not 
mainly emitted in the forward direction. The lower their kinetic energy the 
more isotropic their angle distribution becomes |36| . 

Evidently, there are more interaction processes like ejection of Auger electrons 
or characteristic X-ray photons. But they are very unlikely in the energy range 
considered. 

Although inelastic collisions are decisive for the energy transfer, the radiation 
damage in the patient strongly depends on the spatial distribution of electrons in 
their passage through matter. They dominantly undergo multiple scattering events 
with small deviations. However, single backward scattering events also occur fre- 
quently which leads to tortuous trajectories of electrons. To a big extent such 
trajectories are due to elastic collisions. We therefore focus on very accurate and 
realistic simulation of elastic processes in our model. This is achieved by trans- 
port coefficients extracted from the ICRU 77 database [29]. Inelastic transport 
coefficients are obtained in the same way. 

Elastic and soft inelastic events lead to small energy loss. With kinetic energies 
above 1 keV, electrons are assumed to lose their energy continuously [32]. Because 
of this, we implement the CSD approximation to model energy loss of electrons. 
Hence, we neglect large energy loss fluctuations caused by hard inelastic collisions. 

Bremsstrahlung effects are considered in our model in a very restricted way: 
Only energy transfer of electrons after soft bremsstrahlung collisions are simulated 
by means of the radiative stopping power. However, the effect of hard photon 
emission as well as the transport of photons are disregarded so far. 



2.2. The Generalized Fokker-Planck Approximation for Electrons. Parti- 
cles can be described by a six dimensional phase space (r, £2) G (V^ x / x S*^) 
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with 

spatial variable r — {xi, X2, x^) £ V <Z M.^ open and bounded 

energy variable E el ^ [E^inn, £^max] C M 

direction variable il = (fii, ri2, i^s) e S"^ <Z unit sphere 

= (Vl - /x2cos(0), \/l - /i2sin(0),/i), /i = cos(6'). 

If particle transport takes place in an isotropic and homogeneous medium in which 
interaction processes are Markovian and particles do not interact with themselves, 
their distribution can be described by the unique solution of the time-dependent 
linear BTE. For practical applications in radiotherapy we are faced with issues like 
distributions of the deposited energy in tissue or penetration depths of the beam. 
For many purposes it is therefore sufficient to know the steady solution: 

(1) (Ja-^{r,E,n) + n-V^{r,E,n) = LB^{r,E,n) 
with 'i!{r,E,Q,)^^b{r,E,Q,) W e dV, fl- n < 0, E e I . 

^fci-^jil) is called angular flux and denotes the distribution of particles travel- 
ling in direction f2. Boundary conditions (BC) are imposed on the angular flux 
depending on the incident beam. Absorption of particles is described by the ab- 
sorption cross section a a- The right hand side 

LB'^{r,E,n) := [ [ as{E\ E,n- g!)^{r, E\g!)dn' dE' -i:,{E)^{r, E,n) 

is known as linear Boltzmann Operator and describes scattering. Its integral con- 
tains the differential scattering cross section (DSCS) as{E' ,E,Q^- O!) character- 
ising interaction mechanisms in which particles are deflected. The dot product 
£2 • = cos(0o) — Mo indicates that the scattering probability only depends on the 
scattering angle. This implies that the deflection of scattered particles is axially 
symmetrical around the direction of incidence £2'. Integrating the scattering kernel 
(Js{E' , E,[l-[l') over all angles and energies, one gets the total differential scattering 
cross section (TDSCS) 

(2) ^siE)^2TT / as{E',E,fio)dfiodE'. 

Jo J-l 

The angle and energy integral in the Boltzmann operator is the main difficulty 
in its numerical solution. That is why an important aim is to develop accurate 
approximations. However, up to now there is no predominant method used for all 
types of particles. In fact, depending on specific particle properties, one has to 
choose an appropriate approximation. 

One crucial property of elastic DSCSs in water is a sharp peak in the forward 
direction [301 . To express this mathematically we define the positive n-th scattering 
transport coefficient (STC) 

/•OO pi 

(3) ^n{E):=2nJ J (1 - ^o)"cr,(£;', i;, ^0)^/^0^^^' for all n > 

and assume that as n increases, the coefficients ^,1 fall off sufficiently fast, i.e., 

(4) ^n+i{E) <^ £.n{E) for all n > and £; e /. 

Additionally, elastic scattering often entails small energy loss. A first approximation 
is therefore to expand the scattering kernel around /io = 1 and E = E' . In this 
way Pomraning [35] showed that the already known FP operator is the lowest-order 
asymptotic limit of the integral operator Lb- In [39 this FP operator is derived as 
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a first order angle approximation to Lb'- 
Fokker-Planck Lpp 

(5) where Ls^-® = Xfp^'® + 0{e) for e < 1. 

Since the spherical Laplace-Beltrami operator 



L = 



with /i = cos(6') 



9/i ^ ^ ^ 9/i \ — ^"^ d(j) 

is differential in angle the nonlocal integral Boltzmann operator is now ap- 
proximated by a local differential operator. The crucial point is that an integro- 
differential equation is transformed into a partial differential equation. Although 
discretizations of differential equations often lead to large linear systems their nu- 
merical effort turns out to be much lower. This is due to the local character of 
differential equations which bring along much sparser matrices. 
Pomraning's resulting Fokker-Planck equation for particle transport in an isotropic 
medium yields: 

(6) E.ann. vn,. e, m = ^«(., e, m + '""fc^'/fc^'^" . 

^(r, e) is called stopping power defined by 

(7) S{r,E)= [ [ {E - E')as{E',E,n-n!)dndE'. 

Jo J47r 

The standard Fokker-Planck approximation is a frequently used method to describe 
transport processes in media with highly forward-peaked scattering. Comparisons 
with real data, however, reveal that many scattering processes of interest contain a 
small but sufficient amount of large-angle scattering. To gain higher order asymp- 
totic approximations to Lb one could expand Lb to a 

Polynomial Operator Lp„ := I]m=i a«,mi" with a„^„ e 0{^rn) = 0(e™"^) 

with = ip„«'(fl) + 0{e") for all n>l. 

However, Leakeas and Larsen showed that eigenvalues of Lp„ might become positive 
so that the angular flux could become infinite [SS]. This served as motivation for 
them to develop asymptotically equivalent operators to Lp„ which remain stable 
and preserve certain eigenvalues of Lb- We have summarized the details of the 
computation in the appendix. We end up with Generalized Fokker-Planck (GFP) 
equations which incorporate large-angle scattering and are therefore more accurate 
than the conventional Fokker-Planck equation: 

(8) aa^ir, E,n)+n- V^I/fc E, Q) = LcFpMr, E, Q) + ^iSir, EMl, E,n)) _ 
For positive coefficients ai{E), j3i{E) and m G N, the GFP operators are defined by 
-^GFPa™ := X! oii{E)L{I - (3i{E)Ly^ and Lgfp2„,+i ■= Lgfp2^ + am+i{E)L. 

2.3. Determination of Physical Quantities. The fundamental part in GFP 
theory is based on the assumption of forward-peakedness (eq. Q). Its expected 
accuracy strongly depends on the behaviour of transport coefficients £,n{E) which 
differ for different projectiles and materials. The ICRU Report 77 provides differ- 
ential cross sections for elastic and inelastic scattering of electrons and positrons 
for different materials and energies between 50 eV and 100 MeV [55]. To obtain 
transport coefficients £,n{E) we use these cross sections and proceed in the following 
way: 



APPLICATION OF GENERALIZED FOKKER-PLANCK THEORY 



7 



(a) The ELSEPA code system, distributed with the report, calculates elastic and 
inelastic angular differential cross sections for a fixed energy E 

/>oo 

a^'''"°'(£;,/io) = / af"'°\E',E,fio)dE' 
Jo 

in tabulated form for discrete fiQ and (t'^'''"°'(£', hq). For a predetermined set of 
energy values between 50eV and lOOMeV, data for cr°^''"'=i(i;, ^o) are extracted 
from these files. 

(b) With that we calculate the n-th transport coefficient for a fixed energy E 

^oi,inci(^) = 27rAA - /^o)"^t" ■"°'(^^,A*o)dMo with ^^o = cos(0o) 

via numerical integration of the tabulated cross sections a°^'"^''^{E, /^o) by means 
of the trapezodial rule. Additionally, we multiply the result by the molecular 
density of the transmitted matter 

Na is Avogadro's number, p the mass density of the material and A its molar 
mass. 

(c) Again, all computed results of ''"°'(-^) ^^'^ stored and used as a look-up table. 
To obtain the n-th transport coefficient at the desired energy E this tabulated 
data is linear interpolated. 

Finally, we use the following transport coefficient in our equation: 



Fig. [T| illustrates electron transport coefficients of different order in liquid water. 
Except for the 0-th, every coefficient is strictly monotonically decreasing. For 
E > 10"^, ^1 is always bigger than ^2 but, as E decreases, their deviation re- 
duces more and more. Unfortunately, for increasing n > 2 the difference between 
two consecutive ^„ is so small that the assumption of forward-peakedness is not ful- 
filled. The inelastic transport coefficients are much smaller than elastic transport 
coefficients for n > 1. The higher the order of the transport coefficient the larger 
is the difference between elastic and inelastic ones. 

Our stopping power in eq. Q is equivalent to the physical total stopping power 
as the sum of collision and radiative stopping powers. For different materials both 
are directly included in the files of the ICRU database. Hence, we do a linear 
interpolation to get the stopping power S{E) at the desired energy E. 

3. Deterministic Model for Light Propagation 

3.f. The Generalized Fokker-Planck Equation for Grey Photons. Many 
medical applications like cancer treatment or optical imaging of tumors make use 
of propagation of laser light in tissue. Its behaviour is determined by the solution 
of the steady grey radiative transport equation 



(9) aa^{r,n) + n-v^{r,n)^ 



with -^{r, n) = -^bir, n) yredv,n-n< 0. 



Similar to |2.2| one can derive GFP approximations to the radiative transport equa- 
tion in a straightforward way. The intensity ^(r, £2) describes the radiation power 
flowing in direction CI which is influenced by scattering and absorption coefficients 
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Figure 1: Comparison of scattering transport coefficients in liquid water for 
energies between 50eV and lOOMeV. 



a a and /i^. More important is the scattering kernel CTs(£2 • £2'). It is also charac- 
teristic for biological tissue that this kernel has a sharp peak at • fi' = 1. For 
its simulation mathematically simple scattering kernels with a free parameter are 
used. 



3.2. Models for Scattering Kernels. One often cited scattering kernel is the 
(single) Henyey-Greenstein kernel defined by 

yHG 

(10) af ^(Aio) = -i-fHG{^^o) 

1-^2 

where fHGif^o) wf, ^ r^z^ ^ (-1; 1) 

2(1 - 2.g/io + g-^y/^ 



is the corresponding phase function and Sj^ the TDSCS. The single parameter 
g determines the amount of small- and large-angle scattering. If g ^ 1, erf is 
not only strongly forward-peaked but also includes large-angle scattering. Its value 
depends on the irradiated tissue. For human tissue typical values for g are around 
0.9 (human blood: g = 0.99, human dermis: g = 0.81 [in])- Expanding fnGil^o) in 
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Legendre polynomials gives expressions for ^„ depending only on g: 

(11) a = 1-5 

(12) 6 = ^-25+^3^ 

(13) C3^2-fg + 2g'-lg' 

^ 16 32 32 2 8 3 8 4 

(14) ^4 = y - y g + - ^9^ + 

16 80 200 2 40 3 8 4 8 5 

(15) e5 = y-y5+^3 +7^ "63^ 

This is all we need to calculate ai and /Si for GFP2 to GFP5. It turns out that 
they remain throughout positive. 

To control large-angle and forward-peaked scattering a linear combination of for- 
ward and backward Henyey-Greenstein phase functions was introduced |24j . For 
real constants gi € (—1; 0], g2 € [0; 1), b G [0; 1] and the phase function 

(16) foHGif^o) := bfHG{fJ'0,gi) + (1 - &)/hg(a*o,52), 
the double Henyey-Greenstein scattering kernel is defined by 

yDHG 

(17) ^^(A^o) = ^^/DHG(/i0). 

zvr 

An indicator for the amount of forward or backward scattering is the constant 
b: Setting 6 = the backward scattering phase function fnoilJ-Oi di) vanishes 
whereas b — 1 reduces /d//g(a'o) to the single Henyey-Greenstein phase function 
fHcifJ-Oi 9i)- This provides an opportunity to adapt it to the material of interest. 
Analogous to above procedure for the single HG phase function one can conclude 
that 

(Tsn = 6(5" - 52 ) + 52 

from which all STCs £,n can be calculated. 

4. Numerics for Generalized Fokker-Planck 

4.1. Discretization of Differential GFP Equations. Replacing the right-hand 
side of the Boltzmann equation by the GFP2 approximation operator yields: 

S{r, E) ^^^^'^'-^ + aa^(r, Q) + Q • V^(r, E, Q) = a{E)L{I - (3{E)L)-H{r, E, 0) 

dS{r, E) 



(18) +^{r,E,m- 



dE 



To solve eq. (18) we restate it by setting 

^(°) = ^{r^E,^ and ^^^^ ^ {I ~ 13{E)L)-H'^^\ 
so that it becomes 

5(r, E) ^^^°^Sf;^'-^ -I- aa*(") (r , il) + Q • V*(°) (r , E, Q,) = a{E)L^^^^ (r, E, Q) 
oE 

(19) +^('\r,E,mf-^^ 

(20) {I - l3{E)L)^^^\r,EM) = *^°'(z:, 0)- 

These equations form a coupled system of second-order differential equations with 
the angular momentum operator L. Solving this system requires differencing schemes 
in space and angle. Initial and boundary conditions are imposed on 
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The asymptotic GFP analysis transformed the original BTE into a new type of 
equation which requires an additional condition to the energy variable: 



(21) 



i?oo denotes a large cutoff energy. In the numerical simulations, it should be bigger 
than the energy of all particles from the incoming beam. 

A simplified model to be studied is that of a plate which is infinitely extended in x 
and y directions with a thickness d in z direction. Due to symmetry reasons 

• the angular flux is independent of x and y directions and 

• its direction of motion f2 only depends on 9. 

That is why the initial six dimensional problem has now been reduced to a three 
dimensional one. Although it seems that we only describe a one dimensional object 
in space one should not forget the fact that our slab still remains three dimensional. 
From the mathematical point of view the symmetry of this model, however, leads 
to a one dimensional problem in space which decreases computational costs. 
In slab geometry the aforementioned system reduces to 



(22) 
(23) 



_..,<..(.,i,.,,-e!!Hi£:iO., 

oz 

{I - ^3{E)L,,)¥^\z,E,^JL) = ¥''\z,E,^l), 



where the one dimensional angular momentum operator is defined by 



(24) 



d_ 



We solve this system in two steps: 



(1) Obtain the solution ^'^'^^z, E, ^J.) to eq. ([231. 

(2) Plug *(i)(z,£;,Ai) in eq. (|22| and solve the resulting differential equation. 



To achieve accurate results and lower computation times a high-order scheme was 
implemented [43|: 



(25) 



D 



1 



D 



3 + lf2' 



D 



(1) 



(1) 



J-l/2- 



j+1/2 = Dj-i/2 - "^fJ-jWj with Di/2 = = D 



A/+1/2) 



where fij are abscissas for a Gauss-Legendre quadrature rule with weights Wj . 
With knowledge of the already calculated explicit values ^'P' the right-hand side 



of eq. (22) reduces to a single dependence. The resulting partial differential 
equation is discretized with finite differences in the z-direction. Hence, we end 
up with an ordinary differential equation in the energy variable E. Its solution 
is obtained by the embedded 2nd/3rd order Runge-Kutta MATLAB solver ode23 
solving from the initial condition in eq. (21 1 backward in energy to -E = 0. The 
remaining discretizations are of first order in z and of second order in /i. 
Higher order GFP equations are discretized and solved analogously. However, due 



to more frequent occurrence of and (/ — (3{E)Lfj) the right-hand side of eq. ( 22 ) 
becomes more involved and more systems have to be solved. 
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5. Numerical Results 

5.1. Slab Geometry: HG Kernel. First we neglect absorption and start with a 
simpler form of the GFP equation 

(26) 0- V*(r,il) = iGFP„*(z:,il) 

This is solved in slab geometry with the HG scattering kernel from eq. ( [To| ) for 
selected values of the anisotropy factor g. Symmetry properties mentioned above 
yield the following boundary value problem (exemplarily stated for GFP2): 

(27) /i = aL^^^^^ (z, m) 

(/-/3i^)vI/W(z,M)-*(")(^,/,) 

BC: «'(")(0,Ai) = 10^ -e-io^i^^)' 1>M>0 

For g = 0.8 and g = 0.95 results were computed by time marching with an adaptive 
Runge-Kutta solver until a steady state was reached. Incoming photons moving in 
positive z-direction at z = are simulated by narrow Gaussian peaks around /i = 1. 
Corresponding graphs illustrate steady solutions and use penetration depth in cm 
as X- and J^-^ ^(z, /i, s)dii in l/(Jcm^s) as y-axis. The latter quantity is sometimes 
called energy density and is related to the dose. Discretization parameters in z 
(110 points) and in /i (64 points) direction were chosen large enough to reach 
convergence. Figs. [2]|3] additionally show converged transport solutions generated 
by a discrete ordinates method (DOM) which we use as benchmark in the following. 

Each distribution forms a monotonically increasing function until it reaches a 
maximum and strictly decreases afterwards. There is a large difference between 
FP, GFP2 and the other GFP approximations. However, GFP3-GFP5 are hardly 
distinguishable. For increasing values of g this discrepancy between GFP3 and 
GFP5 becomes bigger whereas results for GFP4 and GFP5 show throughout no 
distance at all. As to the DOM curve we observe that FP and GFP2 give poor 
approximations for small penetration depths. However, GFP3.4^5 give quite good 
solutions throughout the whole penetration range. 

5.2. Single Slab: Light Propagation in Tissue. Gonzalez-Rodriguez and Kim 
studied light propagation in tissue including both forward-peaked and large-angle 
scattering |24j . They examined several approximation methods and especially im- 
plemented GFP2. We want to augment this with results up to GFP5 and compare 
the resulting simulations with the transport solution. We focus on the following 
problem: 

(28) a,M/(z,/x)+M^^^^ = M. 



BC : ^{z = 0, ^i) = e-i°(i-^)' 1 > ^ > 

^{z = d = 2,^i) = -l<fi<0 

It is a slab geometry with a thickness of d — 2mm disregarding any time dependence. 
Its solution enables to compute reflectance i?(/i) and transmittance T(/n) defined 

by 

T{fi) = ^{fi,d) 1>^>0. 
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Penetration z [cm] 

Figure 2: HG kernel in slab geometry with g=0.80. The difference between GFP3-GFP5 
lines is so small that they almost agree here. 

X 10" 
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Penetration z [cm] 

Figure 3: HG kernel in slab geometry with g=0.95. Solutions for GFP3-GFP5 already 
overlap. 
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A. SINGLE HENYEY-GREENSTEIN KERNEL 



In the first run eq. ( 28 1 was solved with discretizations of 64 points in angle and 80 
points in space using GFP approximations with the HG DSCS = crf^'^- Further 
constants were set to 

5 = 0.98 da = O.Olmm^i = 50mm"^ 

REFLECTANCE: Starting at w 0.23 the reflectance slightly increases and 

attains its maximum at /i w —0.5 (fig. |4]). Although in this interval GFP data show 
discrepancies among each other their results are accurate and GFP3 gives the best 
approximation. For /i > —0.5 the transport solution DOM hunches down more than 
GFP functions and hence, the error increases rapidly. Surprisingly, for /i > —0.25 
GFPs-reflectance values are closer to DOM than those of GFP5. Throughout the 
whole interval FP values give a very poor approximation. 

TRANSMITTANCE: It is almost a straight line only bending for small /i. In 
contrast to the reflectance, a more or less constant distance to the transport so- 
lution is always sustained. To the eye, there are no differences between all GFP 
simulations in a wide range. Only for small /i the functions start to deviate and 
GFP3 data give best results whereas FP is inaccurate again. 



B. DOUBLE HENYEY-GREENSTEIN KERNEL 



Taking the amount of large-angle scattering in biological tissue into account Gonza 
lez-Rodrfguez and Kim applied the double Henyey-Greenstein DSCS to simulate 
transmittance and reflectance in liver tissue. The following fit parameters were 
used: 

51 = 0.85 52 = -0.34 b = 0.86. 

51 — 0.85 provides a forward-peak which is not very sharp. In addition, the com- 
bination of 52 — —0.34 and b = 0.86 contains a significant amount of large-angle 
scattering which leads to increasing STCs: 

a = 0.3166 C2 = 0.3916 ^3 = 0.6058 ^4 = 1-0075 ^5 = 1-7388 

In this case our fundamental assumption is not valid which could negatively affect 
our approximations. Moreover, it is important to emphasize that simulations for 
GFP3-GFP5 ran with some negative coefficients a;, (3i. Nevertheless, our code gave 
reasonable results plotted in fig. [5] for discretization parameters of 64 points in ^ 
and 70 points in z direction. 

REFLECTANCE: Fig.|5]shows 'bump head' functions similarly shaped to those 
of the single HG kernel. The x-coordinates of their maxima are, however, shifted 
to the right. Moreover, for large /i different GFP approximations do not match as 
well as they do in fig. |4] In contrast to the single HG kernel GFP3 data give a poor 
approximation whereas GFP5 is the best one among all shown here. Only for sa 
GFP2 is not able to match GFP5. As expected, FP gives even worse results than 
for the single HG kernel. 

TRANSMITTANCE: This time our transport solution DOM is more peaked at 
jjLK,!. Nevertheless GFP gives more accurate results than in fig.|4] A comparison 
between the two best approximations GFP3 and GFP5 yields small differences 
which enlarge near /i « 0. However, the classical FP deviates from our benchmark 
to a big extent. 

Due to the contribution of large angle scattering numerical computations with the 
double HG kernel are more challenging and, in fact, give GFP coefficients which 




Figure 4: Single HG with g=0.98: Reflectance and transmittance of liver 
tissue arising from a slab geometry with thickness d=2mm and conditions in 
eq. (281. Transmittance is plotted in a semilogarithmic scale. 
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Figure 5: Double HG with = 0.85, (72 ~ —0.34: GFP approximations for 
reflectance and transmittance of liver tissue. Transmittance is plotted in a 
semilogarithmic scale. 
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contradict our assumptions. Nevertheless, the GFP results plotted above approxi- 
mate the transport solution very well and are much more precise than those of the 
FP calculations. 

5.3. Slab Geometry: Electron Propagation in Tissue. For dose calculations 
the following GFP equation is to be solved (examplarily stated for GFP2): 

oz oh 
(29) (/ - /^L^)* (i)(z, s) = vl/(0)(z, E, fi) 

BC: vi/(0)(0,ii;,^) = 105.e-20o(i-P)%-50(£;o-£)^ l>^>0,i?e/. 
^'(°)(d,£:,Ai) = -l<fi<0,E'El. 



The initial boundary value problem in eq. ( 29 ) describes the propagation of elec- 
trons through matter with a monoenergetic pencil beam of energy Eq irradiated 
orthogonally to the boundary surface of the material. This is modelled by a product 
of two narrow Gaussian functions around /i = 1 and E — Eq. After computing the 
solution one can calculate the absorbed dose via 

(30) D{r)=— / S{r,E')■^(^\r,^i,E')d^,dE'. 

P[r) Jo J^i 

T is hereby the duration of the irradiation of the patient and p the mass density of 
the irradiated tissue so that D{r) leads to SI unit J/kg or Gy. 
Several test cases were implemented for 5 MeV and 10 MeV beams. As benchmark 
we used solutions of the MC code systems GEANT4 (standard physics package) 
[HIS] and PENELOPE g?]. However, it should be stressed that all physical mod- 
els were obtained independently. The following criteria are generally employed to 
quantify the accuracy of a dose curve |54l : 2%/2mm (pointwise difference within 
2% or 2mm horizontal distance-to-agreement) in homogeneous and 3%/3mm in in- 
homogeneous geometries. 

A. HOMOGENEOUS GEOMETRY 

Characteristic electron dose profiles in a semi-infinite water phantom are shown 
in fig. |6] and fig. [7j First they provide a high surface dose, increase to a maxi- 
mum at a certain depth and drop off with a steep slope afterwards. Solutions for 
GFP4 and GFP5 are omitted because they overlap with GFP3 in our plot. Except 



for GFP2, computations were performed according to eq. (25) (32 points in fi, 350 
points in z) . Due to better results we applied upwind finite difference discretizations 
for GFP2, equidistant in z (400 points) and /i (200 points). All approximations are 
close to each other because GFP transport coefficients ^n{.E) for water do not fall 
off highly enough within our energy interval. All in all, the calculated results agree 
well with PENELOPE and GEANT4. All dose profiles for a 5 MeV beam satisfy 
the 2%/2mm criterion. As we neglect bremsstrahlung the difference to MC com- 
putations becomes bigger for 10 MeV. In fact, the largest FP and GFP2 distance 
to PENELOPE and GEANT4 becomes 3mm at z ~ 5 cm and hence, they do not 
meet the criterion. 

B. INHOMOGENEOUS GEOMETRIES 



Dose calculation is more challenging in parts of the body where materials of strongly 
varying densities meet. Here, large dosimetric differences between experiments and 
predictions exist As deviations of already five percent in the deposited dose 




Figure 7: 10 MeV electron beam: normalized dose in liquid water. 
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may result in a 20% to 30% impact on complication rates [Hj it is of big importance 
to accurately compute the dose in such transition regions. 

Possible clinical applications for electron beams are for example irradiation of 
the chest wall or the vertebral column. To simulate dose curves on the central 
beam we assume that 10 MeV electrons pass three different materials: muscle (0- 
1.5cm), bone (1.5-3cm) and lung (3-9cm). For all results, parameters for Morel's 
discretization |43j were set to 32 points in ^ and 400 points in z. Fig. [s] illustrates 
approximations up to order three because higher order results overlap with the latter 
on that scale. The agreement with the MC dose profile is very satisfactory although 
bigger differences occur for small penetration depths. The dose differences between 
PENELOPE and FP exceed the 3%/3mm hmit only at the boundary z = 0. 

Radiotherapy gains in importance not least because surgical interventions can 
be avoided. Especially sensitive body areas like the brain, coated by cerebral mem- 
branes, are of big interest. Between those membranes there are many voids which 
means that scattering and absorption properties change abruptly. Therefore we 
consider an air cavity irradiated by a 10 MeV electron beam first penetrating water 
(0-4cm), then air (4-6cm) and later water (6-9cm) again. Similar to pure wa- 
ter, GFP2 calculations yield better solutions for equidistant upwind discretizations 
(200 points in /i and 300 points in z). Remaining curves were obtained by Morel's 
scheme (/i-direction: 32 points, z-direction: 350 points). Except for small penetra- 
tion depths all deterministic solutions are very close to each other and demonstrate 
good aggreement with PENELOPE (fig. [9]). Again FP and GFP2 results show the 
best approximations. Larger, but still comparably small, differences between them 
occur in the air region. Except for the boundary value {z ~ 0), the FP and GFP2 
curves fulfill the 3%/3mm criterion. 

6. Conclusions 

We studied practical applications of Generalized Fokker Planck approximations. 
Numerical examples of GFP solutions for the Henyey-Greenstein kernel in slab 
geometry showed more accurate approximations than FP calculations. Further 
test cases for refiectance and transmittance in liver tissue by means of single and 
doulbe HG kernels also revealed GFP3- and GFPs-results closest to the transport 
solution. 

For electron transport, we derived an ab initio model from the ICRU database. This 
model was compared to publicly available MC Codes (PENELOPE and GEANT4) 
which in turn has been benchmarked against experiments. We extracted the stop- 
ping power, elastic and inelastic cross sections from the ICRU database and trans- 
formed them to transport coefficients needed for GFP computations. Dose distribu- 
tions for electron beams were performed without additional coupling to photons and 
positrons. We are aware that our physical model neglects important interactions 
like energy straggling and hard radiative events with emission of photons. They 
are inevitable for accurate dose calculations with high-energy electrons. However, 
in our energy range they are less frequent and regarded as extensions for improved 
models in future. And in fact, comparisons of GFP approximations with Monte 
Carlo calculations reveal dose profiles which agree well in both homogeneous and 
inhomogeneous geometries. 
Several tasks for further examination remain: 

(i) The first step towards real dose calculations from CT data is an extension to 
two space dimensions. As the GFP theory was derived for angular fiuxes in 
3D space this should only be a challenge to the numerical and programming 
approach. 
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Figure 8: Normalized dose curves of 10 MeV electrons irradiated on the back of the body 
(32 points in fi, 300 points in z). 
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Figure 9: 10 MeV electron beam: normalized dose in liquid water with air cavity. 
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(ii) Due to a rising demand for proton therapy facilities the adaption of the GFP 
theory to protons is certainly an interesting subject of further study. 

(iii) To improve computational results it is necessary to include more physical 
phenomena. Especially for high energy electron beams it is inevitable to 
simulate the transport of bremsstrahlung quanta. 
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Appendix B. Derivation of Generalized Fokker-Planck Operators 

GFP2 

If denotes one eigenvalue oi Lb for n > and a, /3 are two positive constants 
the Generalized Fokkcr-Planck (GFP) operator defined by 

LP2 + 0(£') = Lgfp^ ■= aL{I - I3L)-^ 

(35) =aL + aPL'^ + 0{aff) 

will have to satisfy three properties to substitute Lp2 in the favoured way: 

(1) Eigenvalue preservation 

an{n + 1) 



1 + /3n(n + 1) 



\GFP2 L \l 



for 



1,2 



Multiplying above equation by (1 + (3n{n+l)) ^ and dividing by n(n+l) 
we conclude: 
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n(n + 1) 



n= 1,2 



1 -CTal 




a 




'(Jal/2 


1 -<To2 




A 




.(7a2/6 



APPLICATION OF GENERALIZED FOKKER-PLANCK THEORY 



21 



(2) Order ©(a/?^) = 0{e^) 

(3) Equivalence 



a an is a quantity which can be expressed in terms of ^n' 
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In the following item (1) is first transformed to a system of linear equations and 
thereafter solved for the desired GFP coefficients (here: a and /3). However, the 
final equations to be solved are non-linear for GFP operators of order n > 3. In 



this case of order n = 2 eqs. (37l-(38l yield 



a = ^ + ^ and /3 



Going on with item (2) it is to be checked 
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As equivalence condition follows straight forward it has been shown that the oper- 
ator Lgfp2 is an ©(e^) approximation to Lb whose first three eigenvalues agree. 
Now it is quite intuitive to apply a similar procedure to higher order operators. 
We determined explicit solutions for GFP2-GFP5 coefficients and performed ver- 
ifications for items (l)-(3). According to GFP3 all items were checked without 
computer support. The asymptotic behaviour of GFP operators of order four and 
five was, however, checked by means of a symbolic toolbox. To keep our following 
description short we confine ourselves to final results. 



GFP3 



LP3 + O(e^) = Lgfp, ■■=aiL{I - (3iL)~^ + 
(42) =(ai + a2)L + aipiL^ + ai(3fL^ + 0{ail3f) 
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(1) Eigenvalue preservation 
(43) {ai + a2) - aanf^i + n{n+ l)Pia2 

6(27C| + 5C|-2466) 



n(n + 1) 
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(3) Equivalence 



a2 = 



To guarantee that LcfPs = Lp3+0{e^) all coefficients of in eq. (42 1 and 
eq. (32 1 must coincide. Let ai,/3i,a2 be the defining positive coefficients 
of LqfPs as stated in eqs. (44)-(46l then one can show that they satisfy 

I. ai + a2 = — + — + — + ©(£ ) 
ni. a,P^, = ^+0{e'). 

GFP4 

Oie^) = Lgfp, -.^aiLil - p^L)-^ + a2L{I - 132L)-' 
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GFP5 

XP5 + 0(e'') - £gfp5 :-aiL(/ - PiL)-^ + a2L{I - P2L)-^ + a^L 

=(ai + a2 + a3)L + (ai/3i + a2(32)L^ + {ai(3l + a20i)L^ 
(49) +(ai^i3 ^ c^2/3|)i* + (ai/3f + a2/3|)i-^ + 0(ai/3i^) + 0(a2/3|) 

(1) Eigenvalue preservation 
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I. ai + a2 + 0:3 




III. ai/j2 + a2/3| = 

IV. ai/?? + a2/3| = 
V. + a^Pl = 



II. ai/Ji + a2/?2 = 
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All linear and nonlinear equations stated above lead to explicit solutions for ai 
and Equations posed for GFP2 and GFP3 actually deliver uniquely determined 
constants whereas for higher GFP operators there is no guarantee for unique or 
even real valued and Nevertheless, it is important to emphasize that the 
resulting values of ai and (3i must be positive. Otherwise, eigenvalues X^^^'' of a 
GFPfe operator could become negative. For DSCS (7^(0 • 0') of different materials, 
and thus different this has to be checked separately. 
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